{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Damped, driven linear oscillator: various cases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The examples shown in this notebook were taken from M. Tabor, Chaos and Integrability in Nonlinear Dynamics, Chap. 1\n",
    "\n",
    "Here, we consider numerical solutions to ODEs which represent a system subjected to some form of external time-dependent force $F(t)$. Particularly interesting are those cases where $F(t)$ is a periodic function, for example, $F(t) \\approx cos(\\Omega t)$. Consider as an example the damped, driven linear oscillator\n",
    "\n",
    "$$\n",
    "\\ddot{x}+\\lambda \\dot{x}+\\omega^2 x = \\epsilon F(t)\n",
    "$$\n",
    "\n",
    "where $\\epsilon$ can be thought of as a 'coupling parameter' --in the limit $\\epsilon\\rightarrow 0$, the system becomes autonomous again. We can re-write this differential equation as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "        \\dot{x}&=y \\\\\n",
    "        \\dot{y}&=-\\lambda y -\\omega^2 x + \\epsilon \\cos(\\Omega t) \\\\\n",
    "    \\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We will study this system, and consider several distinct cases. We will integrate it using `TaylorIntegration`, and visualize the results using `Plots`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition macroexpand(Module, Any) in module Compat at /Users/Jorge/.julia/v0.6/Compat/src/Compat.jl:1491 overwritten in module MacroTools at /Users/Jorge/.julia/v0.6/MacroTools/src/utils.jl:64.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using TaylorIntegration, Plots\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we start, we fix some integration parameters (the order of Taylor expansions and the local absolute error tolerance):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const abstol = 1e-30\n",
    "const order = 28"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 1. The driven oscillator (non-resonant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will consider the non-resonant, frictionless, linear oscillator:\n",
    "\n",
    "$$\n",
    "\\ddot{x}+\\omega^2x=\\epsilon\\cos(\\Omega t),\n",
    "$$\n",
    "which we will rewrite as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\dot{x}_1 &= x_2 \\\\\n",
    "\\dot{x}_2 &= -\\omega^2x_1+\\epsilon\\cos(x_3) \\\\\n",
    "\\dot{x}_3 &= \\Omega\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This ODE is subject to initial conditions $x(0)=x_0$ and $\\dot{x}(0)=v_0$. Here, we have $\\lambda=0$ (no friction), and $\\omega \\neq \\Omega$ (non-resonant condition)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "drivenosc! (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const Ω = 1.0 #forcing frequency\n",
    "const ϵ = 0.5 #forcing amplitude\n",
    "const ω = 1.1 #\"natural\" frequency\n",
    "const T = 2π/ω #period associated to oscillator's ω\n",
    "\n",
    "#the driven linear oscillator ODE:\n",
    "function drivenosc!(t, x, dx)\n",
    "    dx[1] = x[2]\n",
    "    dx[2] = -(ω^2)*x[1]+ϵ*cos(Ω*t)\n",
    "    nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial time, initial condition and final time are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "571.1986642890532"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const t0 = 0.0 #initial time\n",
    "const x0 = [1.0,0.0] #initial condition\n",
    "#const x0 = [1.0,0.0,Ω*t0] #initial condition\n",
    "const tmax = 100*T # 200*T #final time of integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the particular solution is:\n",
    "\n",
    "$$\n",
    "x_\\mathrm{part}(t)=\\frac{\\epsilon}{\\omega^2-\\Omega^2} \\cos(\\Omega t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x_part (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_part(t) = ϵ*cos(Ω*t)/(ω^2-Ω^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The overall solution is:\n",
    "\n",
    "$$\n",
    "x_\\mathrm{gral}(t) = a\\sin(\\omega t +\\delta)+\\frac{\\epsilon}{\\omega^2-\\Omega^2} \\cos(\\Omega t)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "a &= \\sqrt{ \\left( x_0-\\frac{\\epsilon}{\\omega^2-\\Omega^2}\\right)^2+\\left( \\frac{\\dot{x}_0}{\\omega} \\right)^2 } \\\\\n",
    "\\delta &= \\arctan\\left( \\frac{ x_0-\\frac{\\epsilon}{\\omega^2-\\Omega^2} }{ \\dot{x}_0/\\omega } \\right)\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.3809523809523787, -1.5707963267948966)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const a = sqrt( (x0[1]-ϵ/(ω^2-Ω^2))^2+(x0[2]^2/ω^2) ) #amplitude\n",
    "const δ = atan((ω*(x0[1]-ϵ/(ω^2-Ω^2))/x0[2])) #phase shift\n",
    "\n",
    "x_gral(t) = a*sin(ω*t+δ)+x_part(t)\n",
    "\n",
    "a, δ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check if our expression for `x_gral ` is consistent, using `TaylorSeries`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 9.30148402209536e-17 t - 0.3550000000000001 t² - 1.8757992777892314e-17 t³ + 𝒪(t⁴)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_poly = TaylorSeries.Taylor1([0.0, 1.0], 3) #A Taylor polynomial which represents time\n",
    "x_gral(t_poly) #evaluate x_gral at t_poly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the numerical value of the derivatives of `x_gral` up to any order we want! In particular, its 0-th and 1st order coefficients are consistent, within machine epsilon, with the initial conditions $x_0=1$, $\\dot{x}_0=0$.\n",
    "\n",
    "Now, let's integrate the ODE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.371958 seconds (612.46 k allocations: 56.206 MiB, 4.59% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time tT, xT = taylorinteg(drivenosc!, x0, t0, tmax, order, abstol, maxsteps=50000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the solution represents the sum of two sinusoidal signals with different frequencies and therefore we can observe the phenomenon of beats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x vs t, the first steps\n",
    "firsti =1 # length(tT2)-20\n",
    "lasti =10 # length(tT2)\n",
    "myrange=firsti:lasti\n",
    "lint = linspace(tT[firsti], tT[lasti], 10*length(myrange))\n",
    "plot(\n",
    "lint/T,\n",
    "x_gral.(lint),\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "label=\"analytical\"\n",
    "#leg=false\n",
    ")\n",
    "plot!(\n",
    "tT[myrange]/T,\n",
    "xT[myrange,1],\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "label=\"numerical\",\n",
    "leg=false\n",
    ")\n",
    "scatter!(\n",
    "tT[myrange]/T,\n",
    "xT[myrange,1],\n",
    "label=\"numerical\",\n",
    "leg=true,\n",
    "ms=3.0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x vs t\n",
    "plot(\n",
    "tT/T,\n",
    "xT[:,1],\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now, how does the numerical solution compare to the analytical solution? Well, we can measure this via the difference between the analytical solution and the numerical solution at each time step:\n",
    "\n",
    "$$\n",
    "\\Delta x(t) = x_\\mathrm{num}(t)-x_\\mathrm{gral}(t)\n",
    "$$\n",
    "\n",
    "We measure this in machine epsilons (`eps()`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(945.5, 0.6220322210267172, 0.6220322210265072)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the absolute error in machine epsilons\n",
    "Δx = (xT[:,1]-x_gral.(tT))/eps()\n",
    "#absolute error, numerical solution, analytic solution, at last time step:\n",
    "Δx[end], xT[end,1], x_gral(tT[end])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of the integration, the error in the numerical solution is less than 700`eps()`! Now, let's plot the absolute error as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "tT/T,\n",
    "Δx, #the absolute error in machine epsilons\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"absolute error (machine epsilons)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2. The driven oscillator: resonance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "When $\\Omega \\rightarrow \\omega$, the solution\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x_\\mathrm{gral}(t) &= a\\sin(\\omega t +\\delta)+\\frac{\\epsilon}{\\omega^2-\\Omega^2} \\cos(\\Omega t)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "given before, breaks down. This condition is an example of resonance. In this case, the ODE we are trying to solve is:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  \\begin{split}\n",
    "    \\dot{x}_1 &= x_2 \\\\\n",
    "    \\dot{x}_2 &= -\\omega^2x_1+\\epsilon\\cos(\\omega t) \\\\\n",
    "  \\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and the general solution (for $x_1$) is:\n",
    "\n",
    "$$\n",
    "x_\\mathrm{gral}(t) = a\\sin(\\omega t +\\delta)+\\frac{\\epsilon t}{2\\omega} \\cos(\\omega t)\n",
    "$$\n",
    "\n",
    "Note that, once we know $x_1$, then $x_2$ can be calculated straightforward; $x_3(t)=\\omega t$ by definition. The coefficients $a$ and $\\delta$ are given in this case by:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "a &= \\sqrt{x_0^2+\\dot{x}_0^2/\\omega^2} \\\\\n",
    "\\delta &= \\arctan \\left( \\frac{\\omega x_0}{\\dot{x}_0} \\right)\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "drivenoscres! (generic function with 1 method)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the resonant driven linear oscillator ODE:\n",
    "function drivenoscres!(t, x, dx)\n",
    "    dx[1] = x[2]\n",
    "    dx[2] = -(ω^2)*x[1]+ϵ*cos(ω*t)\n",
    "    nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x_gral_res (generic function with 1 method)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const a′ = sqrt( x0[1]^2+(x0[2]^2/ω^2) ) #amplitude\n",
    "const δ′ = atan((ω*x0[1]/x0[2])) #phase shift\n",
    "\n",
    "x_gral_res(t) = a′*sin(ω*t+δ′)+ϵ*t*sin(ω*t)/(2ω) #overall solution with resonance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to test the consistency of our expression for `x_gral`, we can do so by evaluating it using `TaylorSeries`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 1.0 + 6.735557395310444e-17 t - 0.3550000000000001 t² - 1.3583374080542732e-17 t³ + 𝒪(t⁴)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_gral_res(TaylorSeries.Taylor1([0.0, 1.0], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We're ready to perform a Taylor integration for the resonant driven linear oscillator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "571.1986642890532"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const tmax = 100*T # 200*T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.127499 seconds (518.14 k allocations: 54.571 MiB, 22.92% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time tT2, xT2 = taylorinteg(drivenoscres!, x0, t0, tmax, order, abstol, maxsteps=50000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x vs t, the first steps\n",
    "firsti2 =1 # length(tT2)-20\n",
    "lasti2 =10 # length(tT2)\n",
    "myrange2=firsti2:lasti2\n",
    "lint2 = linspace(tT2[firsti2], tT2[lasti2], 10*length(myrange2))\n",
    "plot(\n",
    "lint2/T,\n",
    "x_gral_res.(lint2),\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")\n",
    "plot!(\n",
    "tT2[myrange2]/T,\n",
    "xT2[myrange2,1],\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false,\n",
    ")\n",
    "scatter!(\n",
    "tT2[myrange2]/T,\n",
    "xT2[myrange2,1],\n",
    "leg=false,\n",
    "ms=3.0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x vs t\n",
    "plot(\n",
    "tT2/T,\n",
    "xT2[:,1],\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")\n",
    "plot!(\n",
    "tT2/T,\n",
    "x_gral_res.(tT2),\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution grows unbounded with time! Physically, this means that the external source is imparting energy to the system, and there are no losses of energy due to any damping ($\\lambda = 0$).\n",
    "\n",
    "Once again we ask: how does this numerical solution compare to the analytical solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(237767.0, 1.000000000053305, 1.00000000000051)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the absolute error in machine epsilons\n",
    "Δx2 = (xT2[:,1]-x_gral_res.(tT2))/eps();\n",
    "Δx2[end], xT2[end,1], x_gral_res(tT2[end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "tT2/T,\n",
    "Δx2, #the absolute error in machine epsilons\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"absolute error (machine epsilons)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, try doing longer integrations for the system above. How does the absolute error behave for longer times of integration? Does it grow systematically? Is there a way to know if the error in our integrations is dominated by roundoff errors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. The damped driven oscillator ($\\lambda \\neq 0$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We return to the original problem:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "        \\dot{x}_1&=x_2 \\\\\n",
    "        \\dot{x}_2&=-\\lambda x_2 -\\omega^2 x_1 + \\epsilon \\cos(\\Omega t) \\\\\n",
    "    \\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "i.e., here $\\lambda \\neq 0$. In this case, the overall solution may be written as:\n",
    "\n",
    "$$\n",
    "x(t) = Ae^{-\\lambda t/2}\\cos(\\nu t + \\Delta)+\\frac{\\epsilon\\cos(\\Omega t+\\Delta_\\mathrm{ext})}{\\sqrt{(\\omega^2-\\Omega^2)^2+\\lambda^2\\Omega^2}}\n",
    "$$\n",
    "\n",
    "where the intrinsic frequency $\\nu$ is\n",
    "\n",
    "$$\n",
    "\\nu = \\frac{1}{2}\\sqrt{4\\omega^2-\\lambda^2}\n",
    "$$\n",
    "\n",
    "and the \"external\" phase $\\Delta_\\mathrm{ext}$ is given by\n",
    "\n",
    "$$\n",
    "\\Delta_\\mathrm{ext} = \\arctan \\left( \\frac{ -\\lambda\\Omega }{ (\\omega^2-\\Omega^2) } \\right)\n",
    "$$\n",
    "\n",
    "We actually have three sub-cases,  depending on the value of $4\\omega^2-\\lambda$:\n",
    "\n",
    "+ $4\\omega^2-\\lambda>0$: sub-damping\n",
    "+ $4\\omega^2-\\lambda=0$: critical damping\n",
    "+ $4\\omega^2-\\lambda<0$: super-damping\n",
    "\n",
    "Here, we will analyze only the first case, but we only need to change the values of $\\omega$ and $\\lambda$ and then run the integration again to analyze those other cases too!\n",
    "\n",
    "The amplitud $A$ and phase $\\Delta$ of the homogeneous part are determined by the initial conditions and are given by\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "A &= \\sqrt{(A\\cos\\Delta)^2+(A\\sin\\Delta)^2} \\\\\n",
    "\\Delta &= \\arctan\\left( \\frac{A\\sin\\Delta}{A\\cos\\Delta} \\right)\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "A\\sin\\Delta &= \\frac{1}{\\nu}( \\frac{\\lambda}{2}(\\alpha-x_0)-(\\dot{x}_0+\\beta) ) \\\\\n",
    "A\\cos\\Delta &= x_0-\\alpha\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and, in turn,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\alpha &= \\frac{\\epsilon\\cos(\\Delta_\\mathrm{ext})}{\\sqrt{(\\omega^2-\\Omega^2)^2+\\lambda^2\\Omega^2}} \\\\\n",
    "\\beta &= \\frac{\\epsilon\\Omega\\sin(\\Delta_\\mathrm{ext})}{\\sqrt{(\\omega^2-\\Omega^2)^2+\\lambda^2\\Omega^2}}\n",
    "\\end{split}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: redefining constant Ω\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.2"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const Ω = 2.0 #forcing frequency\n",
    "const ϵ = 0.5 #forcing amplitude\n",
    "const ω = 1.1 #\"natural\" frequency\n",
    "const λ = 0.2 #damping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "drivdamposc! (generic function with 1 method)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the driven, damped linear oscillator ODE:\n",
    "\n",
    "function drivdamposc!(t, x, dx)\n",
    "    dx[1] = x[2]\n",
    "    dx[2] = -λ*x[2]-(ω^2)*x[1]+ϵ*cos(Ω*t)\n",
    "    nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A=5"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: redefining constant x0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".193145415819222\n",
      "Δ=-0.08222027678306153\n",
      "AcosΔ=5.175602019108521\n",
      "AsinΔ=-0.42650093744798\n",
      "Δext=3.2839914642983628\n"
     ]
    }
   ],
   "source": [
    "#two-variable versions of atan(y/x)\n",
    "myatan(x, y) = y>=zero(x)?( x>=zero(x)?atan(y/x):(atan(y/x)+pi) ):( x>=zero(x)?(atan(y/x)+2pi):(atan(y/x)+pi) )\n",
    "myatan2(x, y) = y>=zero(x)?( x>=zero(x)?atan(y/x):(atan(y/x)-pi) ):( x>=zero(x)?(atan(y/x)):(atan(y/x)+pi) )\n",
    "\n",
    "const t0 = 0.0\n",
    "const x0 = [5.0,0.0] #initial condition\n",
    "\n",
    "const ν = 0.5*sqrt(4(ω^2)-λ^2) #the intrinsic frequency\n",
    "const Text = 2π/Ω #period associated to external frequency Ω\n",
    "const Δext = myatan2( (ω^2-Ω^2), (-λ*Ω) ) #the \"external\" phase\n",
    "\n",
    "#some auxiliary variables... (see text)\n",
    "const α = ϵ*cos(Δext)/sqrt( (ω^2-Ω^2)^2+(λ^2)*(Ω^2) )\n",
    "const β = ϵ*Ω*sin(Δext)/sqrt( (ω^2-Ω^2)^2+(λ^2)*(Ω^2) )\n",
    "const AcosΔ = x0[1]-α\n",
    "const AsinΔ = ( (λ/2)*(α-x0[1])-(x0[2]+β) )/ν\n",
    "\n",
    "const A = sqrt( AsinΔ^2+AcosΔ^2 ) #the homogeneous amplitude\n",
    "#we have to be careful with the homogeneous phase, when inverting tan:\n",
    "#if atan( AsinΔ./AcosΔ ) < 0\n",
    "#    const Δ = atan( AsinΔ./AcosΔ )+π #homogeneous phase, case 1\n",
    "#else\n",
    "#    const Δ = atan( AsinΔ./AcosΔ ) #homogeneous phase, case 2\n",
    "#end\n",
    "\n",
    "const Δ = myatan2( AcosΔ , AsinΔ )\n",
    "\n",
    "println(\"A=\", A)\n",
    "println(\"Δ=\", Δ)\n",
    "println(\"AcosΔ=\", AcosΔ)\n",
    "println(\"AsinΔ=\", AsinΔ)\n",
    "println(\"Δext=\", Δext)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x_ddo (generic function with 1 method)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the general solution to the damped driven linear oscillator:\n",
    "x_ddo(t) = A*exp(-λ*t/2)*cos(ν*t+Δ)+ϵ*cos(Ω*t+Δext)/sqrt( (ω^2-Ω^2)^2+(λ^2)*(Ω^2) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we check for consistency of our expression between the analytical solution, `x_ddo`, and the given initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 5.0 - 7.632783294297951e-17 t - 2.7750000000000004 t² + 0.18500000000000005 t³ + 𝒪(t⁴)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_ddo(TaylorSeries.Taylor1([0.0, 1.0], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final time of integration is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: redefining constant tmax\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "314.1592653589793"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const tmax = 100*Text"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready to integrate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.170335 seconds (523.14 k allocations: 55.079 MiB, 11.53% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time tT3, xT3 = taylorinteg(drivdamposc!, x0, t0, tmax, order, abstol, maxsteps=50000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the solution $x(t)$ look like during the first few steps? Let's plot it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# x vs t, the first firsti3-lasti3 steps, starting from the firsti3-th step\n",
    "\n",
    "firsti3 = 1 # length(tT3)-200 #1 # length(tT3)-20\n",
    "lasti3 = 600 # length(tT3) #10 # length(tT3)\n",
    "myrange3=firsti3:lasti3\n",
    "lint3 = linspace(tT3[firsti3], tT3[lasti3], 10*length(myrange3))\n",
    "plot(\n",
    "lint3/Text,\n",
    "x_ddo.(lint3),\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false\n",
    ")\n",
    "plot!(\n",
    "tT3[myrange3]/Text,\n",
    "xT3[myrange3,1],\n",
    "xaxis=\"time, t\",\n",
    "yaxis=\"x(t)\",\n",
    "title=\"x vs t\",\n",
    "leg=false,\n",
    ")\n",
    "scatter!(\n",
    "tT3[myrange3]/Text,\n",
    "xT3[myrange3,1],\n",
    "leg=false,\n",
    "ms=3.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "As $t$ increases, the homogeneous (exponential) part of the solution decays and one is finally left with\n",
    "\n",
    "$$\n",
    "\\lim_{t\\to\\infty} x(t) = \\frac{\\epsilon\\cos(\\Omega t+\\Delta_\\mathrm{ext})}{\\sqrt{(\\omega^2-\\Omega^2)^2+\\lambda^2\\Omega^2}}\n",
    "$$\n",
    "\n",
    "Note that, unlike the case $\\lambda=0$, although the amplitude of the motion becomes large at the resonance $\\Omega\\to\\omega$, it does not diverge. Physically, this asymptotic solution corresponds to steady oscillations of constant energy, which represents a balance between the energy pumped to the system and the energy dissipated by friction.\n",
    "\n",
    "As always, we ask, just how good, quantitatively, is the numerical solution compared to the analytical solution? To answer that (again, as always), we will plot the absolute error as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(63.75, -0.17560201910850107, -0.17560201910851522)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#the absolute error in machine epsilons\n",
    "Δx3 = (xT3[:,1]-x_ddo.(tT3))/eps();\n",
    "Δx3[end], xT3[end,1], x_ddo(tT3[end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "tT3/Text,\n",
    "Δx3, #the absolute error in machine epsilons\n",
    "xaxis=\"t (natural periods)\",\n",
    "yaxis=\"absolute error (machine epsilons)\",\n",
    "title=\"Absolute error vs time\",\n",
    "leg=false\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0-rc3",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
